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Abstract 

We study the problem of finding probability densities that match given European call option 
prices. To allow prior information about such a density to be taken into account, we generalise the 
algorithm presented in [27] to find the maximum entropy density of an asset price to the relative 
entropy case. This is applied to study the impact the choice of prior density has in two market 
scenarios. In the first scenario, call option prices are prescribed at only a small number of strikes, 
and we see that the choice of prior, or indeed its omission, yields notably different densities. The 
second scenario is given by CBOE option price data for S & P 500 index options at a large number 
of strikes. Prior information is now considered to be given by calibrated Heston, Schoebel-Zhu or 
Variance Gamma models. We find that the resulting digital option prices are essentially the same as 
those given by the (non-relative) Buchen-Kelly density itself. In other words, in a sufficiently liquid 
market the influence of the prior density seems to vanish almost completely. 
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1 Introduction 

Many financial derivatives are valued by calculating their expected payoff under the risk-neutral mea- 
sure. In case this measure is given by a risk-neutral density, then for path-independent derivatives the 
expectation can be obtained by integrating the product of the payoff function and the density. 

If a financial model has been chosen for the pricing in a market in which many derivative products are 
actively quoted, then often, due to the limited number of model parameters, this model will be unable to 
perfectly match the market quotes, and a compromise must be made during model calibration by using 
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some kind of "best-fit" criterion. Of course, having the calibrated model brings the advantage of also 
being able to price path-dependent derivatives. 

If no model has been chosen, one can try to imply from the market data for a given maturity a 
probability density function that leads, by integration as described above, exactly back to the quoted 
prices. However, unless the market is perfectly liquid, there will be infinitely many densities that match 
the price quotes, and some criterion for the selection of the density will have to be applied. One such 
criterion is to choose the density that maximises uncertainty or, in another word, entropy. The idea is 
that between two densities matching the constraints imposed by the market prices, the one that is more 
uncertain - where "uncertain" means, very roughly speaking, spreading probability over a large interval 
instead of assigning it to just a few points, where possible - should be chosen. In general, applying the 
criterion of entropy delivers convincing results. For example, on the unit interval [0, 1], if no constraints 
are given, the density with the greatest entropy - the maximum entropy density (MED) - is the uniform 
density. On the positive real numbers [0, oo[, if the mean is given as the only constraint, the entropy 
maximiser is the exponential density. There is no entropy maximiser over the real numbers K, but if 
the mean and variance are imposed as constraints, then the density with largest entropy is a Gaussian 
normal density. 

The concept of entropy has its origins in the works of Boltzmann [Fj in Statistical Mechanics and 
Shannon [29j in Information Theory, and an important recent application has been by Villani |31) and 
others in the field of Optimal Transport. In Finance, too, entropy has become a popular tool, as a survey 
of recent literature (see for example [3] [7J [8j [9] [TT1 [18j 120]) confirms. 

A third approach, which we will use to combine the two approaches described above, is to take a 
density p, which may be near to matching some imposed constraints, as a prior density and to then find 
a density q that is "as close" as possible to p and exactly matches these constraints. The criterion we 
will employ to measure the "distance" between the densities q and p is that of relative entropy. Since we 
are trying to depart as little as possible from the prior p, our goal will be to find the minimum relative 
entropy density (MRED) q that matches the given constraints. Of course, if the prior p already matches 
the constraints, then we can take as our solution q — p, since the relative entropy of p with respect to 
itself is zero. Relative entropy was introduced by Kullback and Leibler [35] and is also known as the 
Kullback-Leibler information number I = I(q\\p) or I -Divergence. Although it is always non-negative and 
can be used as measure of distance, it is important to stress that it is not a metric in the mathematical 
sense, since usually I(q\\p) =/= I(p\\q), and the triangle inequality is not satisfied. 

In the study we carry out in this paper, the prior density function of the asset price for a fixed 
maturity is given by a model, such as the Black-Scholes model or the Heston stochastic volatility model. 
Depending on the model in question, this prior density will be either directly available in analytical form 
(in the case of the Black-Scholes model a log- normal distribution), or have to be obtained numerically (in 
case of the Heston model via Fourier inversion). The main impact of this will be on computation time, 
but otherwise the difference is of minor consequence. The algorithms we propose to calculate the MRED 
q (with respect to p) satisfying some constraints given by European option prices are slight extensions 
of the two algorithms presented in [25] and [27] . 

In the first case, the option data consists of call and digital call prices (see section[2j, and in the second 
case only of call prices (section [3]) . The presence of the digital prices considerably simplifies the problem: 
if only the call prices are imposed, say n of them, the problem consists in finding the minimum of a real- 
valued, convex function (the relative entropy function) in n variables. If one additionally imposes the n 
prices of digital options at the same strikes, the problem simplifies to a sequence of n one-dimensional 
root-finding problems. The second algorithm makes use of the first one by fixing the set of call prices, 
defining a parameter space f2 for arbitrage-free digital prices, and then finding the unique density in this 
family with the smallest relative entropy. 
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The models we take to generate our prior densities are presented in section^ together with a review of 
the characteristic function pricing approach and corresponding Fourier transform techniques. In addition 
to the two models already mentioned above, we also consider the Schoebel-Zhu stochastic volatility model 
and the Variance Gamma model. 

In section [5] we study two market scenarios. In the first one, we take a log- normal Black-Scholes and 
a Heston density as our priors and calculate the MREDs that match given call prices. We then compare 
option prices to those obtained with an MED, and also to those obtained with another log- normal density 
that matches the constraints, and observe that the price differences can be substantial. In the second 
scenario, we take S & P 500 call option prices from the CBOE. This market is very liquid, and for the 
maturity we consider we have quotes for a large set of strikes. We calibrate a Heston, Schoebel-Zhu and 
Variance Gamma model to this data and use the densities generated by these models as prior densities. 
Then, we calculate the three MREDs for these priors, and compare the digital option prices they give 
to those given by the original models, those given by an MED, and finally the market prices themselves, 
which are available in this case. We observe that it makes almost no difference which model is chosen 
for the prior, and that all three MREDs essentially agree with the MED. Finally, section [6] contains the 
conclusion of our article. 



2 Relative Entropy and Densities Matching Call and Digital 
Option Prices 



We begin this section by reviewing the concept of relative entropy, which can be regarded as a way of 
measuring the "distance" between two given densities. Our goal is to apply this measure to the following 
problem: Given a prior density p, coming for example from a model that fits well, but not exactly, 
European option prices observed in the market, how can we deform this density in such a way that it 
exactly matches these prices, but stays as close as possible to the original density under the criterion 
of relative entropy. The options considered here are necessarily European, so that the density is indeed 
sufficient to price them by calculating the discounted expected payoff via integration. 

Next, we show how the algorithm presented in [25] is efficiently implemented in case a prior density p 
is given and the observed prices are those of call and digital options for the same maturity T and the same 
strikes K±, ...,K n . We do not assume that p is given analytically, and therefore the implementation re- 
quires numerical integration. However, the availability of the digital prices allows for an efficient solution 
locally in each "bucket", i.e., interval [Ki,Ki+i[, via a one-dimensional Newton- Raphson rootfinder. 

Finally, we state the results in case the prior is given by Lebesgue- measure, i.e., p = 1, and show how 
all integrals given previously reduce to simple analytic expressions. 

2.1 Relative Entropy 

For two probability distributions Q and P the relative entropy (also called the Kullback-Leibler infor- 
mation number [32] or I-Divergence [13]) of Q with respect to P is defined by 




(1) 



where dQ/dP is the Radon-Nikodym derivative. 
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From the inequality x In x > x — 1 we have 

H(Q\\P)=J§^dP> J {^-l)dP = jdQ-jdP = V VQ. 

We also have i?(Q||P) = if, and only if, Q = P. However, relative entropy is not a metric, since, 
in general, H(Q\\P) ^ H(P\\Q), and the triangle inequality is not satisfied cither. Even the symmetric 
function H (Q\\P) + H (P\\Q) does not define a metric, since it still does not satisfy the triangle inequality 

The Csiszar-Kullback inequality [H |T3] relates relative entropy H to distance between densities in 
the sense of the L 1 (0,oo) norm: 

/•OO 

\\q-p\\ L i:= \q(x)-p(x)\dx<^2H(q\\p), 
Jo 

where 

is the same definition as ([T]) above in terms of densities, which means in particular that convergence in 
the sense of relative entropy implies L 1 -convergence. 



2.2 Finding the Density Matching the Option Data for a Given Prior Density 

We now give a precise formulation of the minimisation problem that we want to solve. Let p be the prior 
density on [0, oof. We assume that p > almost everywhere in [0, oof. We are given a fixed maturity T, 
strictly increasing strikes Kq = 0, K\, K n , K n+ \ — oo, and undiscounted call and digital prices 

Ci := C(K h T)/DF(0,T), A := D{K U T)/DF{Q,T) 

at these strikes. The payoffs of the call and digital options are given in equations ([3]) and DF(0,T) 
denotes the discount factor. Throughout we make the convention C n +i = Ai+i = K n+ iD n+ i = 0. 

Assuming risk neutral pricing, we will determine a density q for the underlying asset price at maturity 
S(T) which minimises relative entropy H{q\\p) under the constraints 

/•OO 

W [(S(T)- Ki)+] =C t , i.e., / {x - Ki)q{x)dx = Ci (3) 

and 

y* oo 

E 9 [l{S(T)> Ki }] = Di, i.e., / q{x)dx = A (4) 

for all i — 0,...,n. In particular, these two constraints for i = mean that q is a density, since 
J q(x)dx = Dq = 1, and that the martingale condition 

/•OO 

E 9 [S(T)] = / xq(x)dx = C 
Jo 

is satisfied, since Co is the forward price of S for time T. 

Formally applying the Lagrange multipliers technique, it can be shown as in [9] that if q minimises 
relative entropy in respect to p, then the Radon- Nikodym derivative g := dQ/dP = q/p is piecewise 
exponential. Therefore on each interval [Ki, the density q is given by 

q(x) = g(x)p(x) = a l e PiX p{x), (5) 
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where on,Pn 6 K, a, > are parameters that still have to be determined using the following two con- 
straints imposed by the option data, which are an equivalent reformulation of the constraints © and 
(|3]) given above, but allow for an easy solution. 

The first constraint follows directly from ^ and is given by 

e^' x p(x)dx = A — A+i, 



from which we have 

A - A+i 



Ik* +1 e ^ iX p( x )dx 



(6) 



'A", 

The second constraint also follows directly from and (U]) and is given by 



pK z . 

a, / xe^ x p(x)dx = (Q + If* A) - (Ci+i + A+i), 
Jic, 

which is the undiscounted price of an "asset-or-nothing" derivative that pays the asset price itself if it 
finishes between the two strikes Ki and Ki + \ at maturity and zero otherwise. Substituting ui from ^ 
gives 

IkT ^ x p(x)dx = (Q + ^A) - + K l+1 D l+1 ) 

jKi+x e ^x p ( x)dx A - A+i 

which we use as an implicit equation for fa. 

We shall show later fCorollarv l2.5[) that, under non-arbitrage conditions, such /3j exists but, firstly, 
we rigorously show that g, given by ([5]) , (|6j) and (0, is indeed a relative entropy minimiser with respect 
to the prior density p. 

Theorem 2.1 For all i = 0, . . . , n, let en and fii be defined by equations ^ and ([7]). Then q : [0, oo[— > R 

given by 

q{x) = a t e^ x p(x) Vx G [ 

minimises H{q\\p). 

Proof This is shown just as Theorem 2.6 in [26 by using Theorem 2.5 by Csiszar also stated there. □ 
We define the cumulant generating functions Cj by 

*{P) := In (^J* Z+1 e (Sx p(x)dx^j . (8) 

Proposition 2.2 For i < n, a is twice differentiable and strictly convex. Moreover, its first and second 
derivatives are given by 



Jk' +1 xe l3x p(x)dx 
Ik' +1 ^ x p(x)dx 



4(P) = -rkz, „ ' ( 9 ) 



f5 I+1 x 2 e^p(x)dx f* i+1 eP*p(x)dx - ( f5 l+1 xe^p{x)dx 
c'HP) = — — ±- '-. (10) 

(Ik! +1 e px p(x)dx 
Additionally, if c„(/3) is finite, then the same holds for i = n. 



5 



Proof We can directly apply Theorem VII. 5.1 (p. 239 in [H]), which says that the cumulant generating 
function (called the free energy function there), under the assumption that it is finite, is a real analytic, 
closed, convex function. 

In our case, the function Cj, for i < n, is finite, since p € L x (Ki, ifi+i) and the exponential function 
e 0x belongs to L°°(K,, K i+1 ). 

The theorem says further that derivatives of all orders exist and can be calculated by differentiating 
under the integral sign, which is how we obtain ([§]) and (flU)) . □ 

Notice that in the last and unbounded bucket [K n ,K n+ \[ — \K n ,oo\ the parameter /? is allowed to 
be positive, but we must assume that the product e /3x p(x) entails a finite integral. Notice also that if 
e@ x p{x) belongs to L x {K n , oo), then so does e@ x p(x) for (3 < (3. 

Introducing 

R , = (Cj + KjDj) - (C i+1 + K i+1 D i+1 ) 

A-A+i 

and using the first derivative of the function Ci, we can rewrite ((7]) in the simpler form 

c'm = Ki. (12) 



This last equation is easily solved for with the Newton-Raphson method using (|9J) and f(T0[) . Once 
the density q has been obtained in this manner, i.e., on, (3i have been calculated for i = 0, . . . , n, we can 
calculate European option prices using numerical integration. 

The next results give the existence and uniqueness of such /3j and, consequently, a*. 

Lemma 2.3 Let m E {0, 1}. Then, for any K £ }Ki, we have 

J* x m e^p{x)dx f^ +1 x m e^p{x)dx 
lim — jr- 2 = lim „ = 0. 

J K i+1 x m eP x p(x)dx J Ki x m el 3x p(x)dx 

Proof We shall consider only the limit when j3 — > oo since the other is treated analogously. Choose L 
and M such that K < L < M < Ki + \. Then, we have 

q Ik x m e l3x p{x)dx J K x m e l3x p{x)dx 
fx' +1 x m e^ x p{x)dx J^ 1 x m e l3x p(x)dx 

Applying the First Mean Value Theorem for Integration yields x\ € [Ki, K] and X2 € [L, M] such that 

j-K pK pM pM 

e^ 1 / x m p(x)dx= / x m e f)x p{x)dx and e^ 2 / x m p(x)dx = / x m e f}x p{x)dx. 

JKi JKi JL JL 

Therefore 

< /* x m e^ x p{x)dx ^ /* x m p{x)dx = 
~ Ik +1 x m eP x p{x)dx ~ e@ X2 J™ x m p{x)dx 
where C > does not depend on j3. Since x\ — X2 < K — L < 0, the result follows. □ 

Proposition 2.4 We have lim^_ i . 00 c^(/3) = Ki+i and lim^-y-oo c^(/3) = Ki. 
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Proof Here again, we only consider the first limit since the other is treated analogously. For all 
K € ]Ki, Ki+i[ we have 



4W)= ( , v 



J* * +1 xe l3x p(x)dx 
Ik' +1 el3x p( x )dx 



(/£ xe^p(x)dx S j (f^ 1+1 xe^p(x)dx S j * + 
/* e0*p(z)da;) e^p(x)cte) 1 + 1 



Ik' +1 el3x p( x )dx 



Using the previous Lemma, we obtain that the term inside square brackets goes to 1 as j3 — > oo. Now 
we shall consider the last term above and show that, by choosing a suitable A, it is as close to -K^+i as 
we want. 

Firstly, we assume i < n. Then < oo and, given a small e > 0, we choose if = Ai+i — e. The 

First Mean Value Theorem for Integration gives x\ € [-Ki+i — e, ifi+i] such that 

"Ki+i a„ , \ , rKi. 



I K !+i-e XePX p( X ) dx _ X lJ K !+i-e el3X p( X ) dx 



fxi+x-e ^ x p{x)dx IkIH-c e px p(x)dx 
Now, for i = n, we have -fQ+i = oo and 

fZ° xe@ x p(x)dx A J™ e^ x p(x)dx 



> 



Ik e ^ x p( x )dx f^ e@ x p(x)dx 
Hence, the term above goes to ifj+i = oo as A goes to oo. 



A. 



□ 



Corollary 2.5 For all i = 0, . . . , n, under the non- arbitrage condition Ki < Ki < A.; + i, where Ki is 
defined in there exists a unique Pi £ K such that c^(/3j) = ATj. 

Proof Proposition 12 . 21 gives that c[ is continuous and the last Proposition states that lim^-oo c^(/3) = 
Ki and lim^oo c'i(f3) = Ki+i. Hence the existence follows from the Intermediate Value Theorem. 
Additionally, Proposition 12.21 also gives that is strictly increasing and the uniqueness follows. □ 



2.3 The Special Case of a Lebesgue Prior 



In case the prior p is given by Lebesgue measure, i.e., p = 1, equation ([5J reduces to the following analytic 
expression: 

> e /3K i+1 _ e f3Ki 



In 



for i < n and /? ^ 0, 



\n(K i+ i — Ki) for i < n and /? = 0, 



In - 



and the first and second derivatives ([9]) and (fTUf reduce to 

f A; +1 e^+* - K ie P K * _ 1 

~ /3 



4(0) = 



+ A 4 



A', 



1 

j8 



for i = n and /3 < 0, 

for i < n and /3 7^ 0, 
for i < n and /3 = 0, 
for i = n and j3 < 0, 
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and 



-{K i+ i - Ki) tnBKi ^ „sk;\2 + ln for i < n and f3 ^ 0, 

for i < n and j3 = 0, 



(e0tfi+i - e ^) 2 /3 2 
(K l+1 - Ki) 2 



12 



for i = n and /3 < 0. 



1 

Using these expressions instead of ([8]). (0), ([TU]) allows for numerical integration to be avoided in an 
implementation in this case. 



3 Relative Entropy and Densities Matching Call Option Prices 

Buchen and Kelly describe in [S] a multi-dimensional Newton-Raphson algorithm to find the maximum 
entropy distribution (MED) if only call prices are given as constraints. The entropy H of a probability 
density q is given by 

H(q) = -jq( x ) In ,(*)<*, 

The minus sign in the definition ensures that H is always positive for discrete densities (where the 
integral sign is replaced by a sum in the definition). For continuous densities, H is usually, but not 
always, positive. For example, the uniform density q(x) = u over the interval [0, i] has negative entropy 

H (q) = — J^ u u In u dx — — In u < for u > 1 . 

In [57] , we show how the results of [53] together with the Legendre transform can be applied to obtain 
a fast and more robust Newton-Raphson algorithm to calculate the Buchcn-Kclly MED. The main reason 
the algorithm is more stable is that the Hessian matrix has a very simple tridiagonal form. In section 
II. A of their paper, Buchen and Kelly also consider the case of "minimum cross entropy" (which we call 
relative entropy here) for a given prior density. We now show how essentially the same algorithm as that 
in [37J can be applied to the relative entropy case. We generalise the results of section 4 in [37J describing 
the Entropy H, and the gradient vector and Hessian matrix to the case in which a prior density p is 
given. 

Arbitrage free digital prices must lie between left and right call-spread prices, i.e., 

Ci — Ci-i ~ Ci + i — Ci 
-K^K- 1 >Di> -K^K i > V. = l,...,n. (13) 

We introduce the set Q C R™ of all D = (L>i,...,D n ) e R™ verifying ([l3|). Note that fl is an open 
n-dimensional rectangle. Define qf) as the density obtained as in Theorem 12.11 for given (undiscountcd) 
digital prices D. 

Proposition 3.1 For all D € f2 the relative entropy of q^ with respect to p can be expressed as 

n n 

H{qo\\p) - J2 PilnPi + J2 PiC i^> 

1=0 4=0 

where c* is the Legendre transform of a, Pi := Dj — Di + \ and Ki is given by (jlip . Notice that pi and 
Ki are given purely in terms of option prices, for all i — - 0, n, and so is H(qf ) \\p). 
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Proof Let g^ := q D /p be the piecewise-exponential Radon- Nikodym derivative given in Theorem 12. II 
We have 

H {qb\\P)= I gp(z)m D dx = lnaj / q £) {x)dx + Pi I xq D (x)dx, 

J Ki P\ x ) JKi J Ki 

since g D (x) = ctie@ iX on [Ki, Ki + i[. Then using the relationships 

a i =p i e- c ^K (14) 
d i (Pi) = K ii (15) 

where ([T4"]) is a reformulation of © using ([5]), the proof goes through exactly as the proof of Theorem 
4.1 in [27]. □ 

Note that if the prior density p already satisfies all of the given constraints, then a, = 1 and (3i = for 
all i = 0, n. It follows from the relationship c*(Ki) — PiKi — Ci(/3i) that in this case c*(Ki) = —Ci(pi). 
Since liapi — lno^ + Cj(ft) = Ci(Pi), Proposition 13.11 savs that 7J(g^,||p) = 0, as expected. 

Proposition 3.2 ^4s a function of digital prices, H : Q — > K is differ entiable and, for all D £ Q, we 
have 

3H 

-~-{D) = \ng b {Ki+) - Ing^Kt-), Mi = 1, n, 
oDi 



where g b := q b jp 



g^{Ki-):= lim g 5 (x) = a^ie^ lKi and g D (Ki+) := lim g f) (x)^a i e l3 ' Ki , 
with ai and Pi given by (|14j) and (|15l) . respectively, for all i = 0, n. 

Proof The proof of Theorem 4.2 from [27] carries over exactly as given there by using the generalised 
versions of Ci and c* introduced above instead of the original Ci and c* , and observing the absence of the 
minus sign in the definition of relative entropy. □ 

From this it follows that, if the prior density p is continuous and D minimises H, then the MRED 
Id = 9dP 1S a -^ so continuous. If p is discontinuous, then q will have the same discontinuities as p. 

Proposition 3.3 The function H : fi — >■ M. is twice differentiable and its Hessian matrix at any Del! 
is symmetric and tridiagonal with entries given by 

obf ' Pi-i Pi + pi-^utfi-i) vi<m ' 



d 2 H - 1 , (Kj-KiHKi+i-Kj) 

-(£>) = 1 — , Mi = 1, ...,n - 1, 



dDidDi+i ft p l c"{P l 

where pi = Di — -Dj+i, and Ki and Pi are given by (|lip and (|12[) . respectively, for all i = 0, n. 

Proof Again, the proof of Lemma 5.1 from [37] carries over without difficulty. □ 

Using these formulations, essentially the same Newton-Raphson algorithm as in [27J can be applied 
to find the relative entropy minimiser q. The only differences are that the functions c" must be replaced 
by their relative entropy versions (fTO]) in the Hessian matrix of Proposition ^. 31 and that in each iteration 
step, for a given set of digital prices, the algorithm of the previous section [5] must be used to calculate 
the MRED, instead of its non-relative version. 
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4 Probability Densities for Characteristic Function Models 



In this section, we look at four models that are popular in equity derivatives pricing. Our aim is to use 
the densities they give for the stock price at a fixed maturity as prior densities. In two of the models 
we have chosen, the Black-Scholes model and the Variance- Gamma model, the density is analytically 
available. In the other two, the Heston and the Schoebel-Zhu stochastic volatility models, it is not. We 
therefore give a brief overview of these models and show how to calculate their densities in each case. 

Let p be the density of x(T) := In S(T). To simplify notation, we will usually write x and S, 
respectively, when it is clear from the context that we have fixed the maturity T. Then the density p of 
S itself is given by 

P(S) = g pQnS), 

since J^^p(x)dx = p(ln S) -^dS = p(S)dS by change-of- variables formula. 
If the characteristic function (c.f.) <f> of p, given by 

/oo 
e mx p(x)dx, (16) 
-oo 

is known, as in the Heston [2T] or Schobel-Zhu [5H] stochastic volatility models, then the density of the 
log-spot price can be obtained via Fourier inversion. Since p is a real-valued function, it follows from 
(|T(3|) that <p(—u) — <f)(u), and we have 

P(x) = ±- f°° e-™<l>{u)du 

27T J-00 

l r°° . i r°° . 

= — I e~ lux (j){u)du+ — / e lux (f)(-u)du 



>(«)] du, (17) 



where 3?[z] = [z + z)/2 denotes the real part of a complex number z. It can immediately be seen that 
an anti-derivative of p is given by 



-j^ /*oo g— iux 

Po(x) = -— I — — <j>{v)du. 
2tt iu 



Furthermore, it can be shown in a similar way as the Fourier Inversion Theorem itself, that lim^-y-oo Po(x) 
— ^ and lim-r^oo Pq{x) = i, and therefore the function 



1 1 /*°° p~^ ux 

P{x) = - - - / —^u)du 

2 2tt J_ 00 iu 

lux (j){u) 



3? 



du 



_ 1 1 

~ 2 ~ 7T 

gives an expression for the distribution function. 

For pricing, we use the general formulation of Bakshi and Madan [JJ. This can be used for a large 
class of c.f. models that contains the Heston, Schoebel-Zhu and VG models (see section 2 in [1], in 
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particular Case 2 on p. 218). We have S > K if, and only if, x > \nK. Let 



III := 1 - P(lniT) = - + - / SR 

2 7T 



oo Y„—iu\nK, 

rfu, (18) 



II 



1 1 f 00 

n 2 :=l-P{\nK) = - + - 9? 

2 7T J 



du, (19) 



represent the probabilities of S finishing in-the-money at time T in case the stock S itself or a risk-free 
bond is used as numeraire, respectively. From (fTB)) we can see that </>(—«) = E [e x ] = E [S], so that the 
quotient 



0(-») y.^ e[5] 

contains the appropriate change of measure. 

The price C of a European call option on a stock paying a dividend yield d is then obtained through 
the formula 

C = e- dJ 'SIli ~ e~ rT KTl 2 , (20) 
and the price D of a European digital call option prices through 

D = e- rT U 2 , (21) 

where r is the risk-free, continuously-compounded interest rate. 



The integrals in f| 18[) . (1191) must of course be truncated at some point a, which depends on the decay 
of the c.f. of the model considered. For Heston and Schoebel-Zhu, we take a = 250, and for Variance 
Gamma a — 1000. 

4.1 The Black-Scholes Model 

Let the parameters r, d and T be given as above, and let a > be the volatility of the stock price S(t). 
In the Black-Scholes model, the logarithm x(t) = \nS(t) follows the SDE 

dx{t) = (r-d- ^<r 2 )dt + adW(t) 

Define /i := In 5(0) + (r — d — ^a 2 )T. The density of x — x(T) is normal and given by 

P(x) = -==e-^r (22) 
V 27rcHi 

The characteristic function of p has a very simple form and is given by 

4>( u ) = e ™M-^ 2 « 2 T • ( 23 ) 

Of course it is faster to use (j2"2"j) directly instead of (|2"5|) and (IT7|) . but comparing these two methods 
lets one measure the additional computational burden. 
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4.2 The Heston Model 



One of the most popular models for derivative pricing in equity and FX markets is the stochastic volatility 
model introduced by Heston J5TJ. Let x(t) := In S(t). The model is given, in the risk-neutral measure, 
by the following two SDE's: 

dx{t) = ( r -d- ^v(t))dt + ^{fydWxit), (24) 

dv(t) = (k9 - (k + \)v(t))dt + a^/v(€)dW 2 (t), (25) 

where (dW\(t), dW^it)) — pdt. The variance rate v follows a Cox-Ingersoll-Ross square- root process [11] . 

The parameter A represents the market price of volatility risk. Since we are only interested in pricing, 
we always set A = in what follows (see Gatheral [TI5], chapter 2). 

Heston calculates the characteristic function solution, but as pointed out in [28], there is a (now 
well-known) issue when taking the complex logarithm. To be clear, we therefore give the formulation of 
the c.f. that we use. 



Define 



b — ipau + d 2 



b := k + A, d\ := \J (ipau — b) 2 + a 2 (f>(i + u), d 2 ■= —d\ and g := — . (26) 

o — ipau — d 2 



Introducing 



k6 ( 1 — ae d2T 
C :={r- d)uiT + — (6 - ipau + d 2 )T + 2 In 



D := 



b — ipau + d 2 1 



a 2 1 - ge^ T ' 

the c.f. of x(T) is then given as 

0(«) =e C+Dv o +iu\nS(0)_ (27) 

Since implementations of the complex square-root usually return the root with non-negative real part 
(di), the key is simply to take the other root (d 2 ), as is done in equation ([26]) . As shown in [I] and [23] . 
this takes care of the whole issue. 



4.3 The Schoebel-Zhu Model 

The Schoebel-Zhu model [5H] is an extension of the Stein and Stein stochastic volatility model [30] with 
correlation p ^ allowed (see also [TU], [33] )• It is described, in the risk- neutral measure, by the following 
two SDE's: 

dx{t) =(r-d- ^v 2 (t))dt + v(t)dWi{i), (28) 

dv(t) = k(6 - v(t))dt + adW 2 {t), (29) 

where (dW\(t) , dW 2 (t)) = pdt. The volatility v follows a mean-reverting Ornstein-Uhlenbeck process. 

The c.f. for this model is given by Schoebel and Zhu in 28 . As Lord and Kahl point out ([23]. 
section 4.2), similar attention has to be paid when taking the complex logarithm in this model's c.f. as 
in Heston's. By directly relating the two characteristic functions (eq. 4.14), they show how the Schoebel- 
Zhu model can also be implemented safely. In the case study in 15.21 presented in the following section 
with SPX option data and a maturity of less than half a year, however, we observed no problems with 
the c.f. originally proposed by Schoebel and Zhu. 
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4.4 The Variance Gamma Model 



The Variance Gamma (VG) process was introduced in [53], [25]. The density for x = In S(T) is given 
explicitly in Theorem 1 in [24] . Define 

lu := — ln(l — Ov — -^ij 2 v) and x — x — In S(0) — (r — d + lu)T. 
Then the density p is given by 



SexpffM . \ Kt ALU(^ + A\, (30) 

I /r/»V5F<rr(5) 1^+0^ V" 2 V V" )) 

where Y is the Gamma-function and K is the modified Bessel function of the second kind. 

If the parameter v is set to zero, the c.f. of p reduces to the Black-Scholes c.f. given in ([23|) . Otherwise, 
it is given by 

_ x. 

4>(u) = e <»ln(S(0))+(r-*fw)T A _ + I^y 2 ^ " . (31) 

Lord and Kahl show ([23], section 4.1) that this formulation of the c.f. is safe. 

Again, as with the Black-Scholes model, since both the density and the c.f. are available, it is possible 
to compare the two different methods (pOl) vs. (I3T1) and (fl~7|) . 



5 Two Numerical Examples 

5.1 A Fictitious Market and Black-Scholes and Heston Prior Densities 

In our first example, we take a hypothetical market with r = d = 0,T = l,S = F = 100, in which call 
option prices are given by the Black-Scholes formula with volatility a = 0.25. As prior densities, we take 

• PBS, a Black-Scholes log-normal density, but this time with volatility a p = 0.20, 

• ph, a Heston density, with parameters k = 1,9 = 0.04, p — — 0.3, a — 0.25, uo = 0.04, which 
leads to implied volatilities of 0.2418,0.2125,0.1923,0.1855,0.1884 at strikes 60,80,100,120,140, 
respectively. 

We calculate the Buchen-Kelly MED (Lebesgue prior), using the algorithm presented in [27], and the 
two Buchen-Kelly MREDs with priors pbs and Ph, using the generalised algorithm presented in section 
[3J and compare the resulting call and digital option prices to see the influence of the priors. 

The different call and digital option prices are reported in table [T] The densities were calculated 
using call prices at strikes 



• K 


= 0,Xi 


= 100 






• K 


= Q,Ki 


= 60, K 2 


= 100, K 3 


= 140 


• K 


= 0,K! 


- 60, K 2 


= 80, K 3 = 


= 100, K 4 = 120, K s = 140 
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Table 1: Comparison of Call and Digital Prices 







Call Prices 






Digital Prices 




Strike 


MED BK 


MRED BS 


MRED Hcston 


MED BK 


MRED BS 


MRED Heston 


20 


80.0538 


80.0000 


80.0000 


0.9936 


1.0000 


1.0000 


40 


60.3244 


60.0000 


60.0094 


0.9766 


1.0000 


0.9979 


60 


41.1698 


40.0637 


40.3043 


0.9316 


0.9841 


0.9595 


80 


23.5389 


21.9716 


22.5717 


0.8124 


0.7758 


0.7828 


100 


9.9476 


9.9476 


9.9476 


0.4962 


0.4420 


0.4763 


120 


3.6684 


3.6071 


3.3294 


0.1830 


0.2039 


0.1977 


140 


1.3528 


1.0596 


1.0051 


0.0675 


0.0693 


0.0593 


160 


0.4989 


0.2688 


0.3239 


0.0249 


0.0192 


0.0175 


180 


0.1840 


0.0621 


0.1171 


0.0092 


0.0047 


0.0056 


20 


80.0000 


80.0000 


80.0000 


1.0000 


1.0000 


1.0000 


40 


60.0015 


60.0003 


60.0012 


0.9997 


0.9998 


0.9996 


60 


40.1454 


40.1454 


40.1454 


0.9669 


0.9753 


0.9715 


80 


22.5812 


22.0890 


22.3433 


0.7743 


0.7818 


0.7770 


100 


9.9476 


9.9476 


9.9476 


0.4646 


0.4424 


0.4633 


120 


3.7041 


3.7051 


3.5189 


0.1945 


0.1976 


0.1926 


140 


1.2139 


1.2139 


1.2139 


0.0705 


0.0707 


0.0617 


160 


0.3800 


0.3569 


0.4669 


0.0221 


0.0227 


0.0207 


180 


0.1190 


0.0961 


0.2067 


0.0069 


0.0065 


0.0077 


20 


80.0001 


80.0000 


80.0000 


1.0000 


1.0000 


1.0000 


40 


60.0033 


60.0002 


60.0014 


0.9994 


0.9999 


0.9996 


60 


40.1454 


40.1454 


40.1454 


0.9726 


0.9727 


0.9726 


80 


22.2656 


22.2656 


22.2656 


0.7794 


0.7781 


0.7804 


100 


9.9476 


9.9476 


9.9476 


0.4510 


0.4499 


0.4510 


120 


3.7059 


3.7059 


3.7059 


0.1971 


0.1961 


0.1958 


140 


1.2139 


1.2139 


1.2139 


0.0700 


0.0711 


0.0689 


160 


0.3834 


0.3545 


0.4105 


0.0221 


0.0227 


0.0211 


180 


0.1211 


0.0948 


0.1564 


0.0070 


0.0064 


0.0071 
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as respective constraints. These strikes are the ones in boldface in table [TJ 

We see that in the first case, where we had only the forward and an at-the-money call option as 
constraints, there are significant differences in both call and digital prices. The presence of the log- 
normal prior density makes MRED BS call prices cheaper compared to the original BS prices. Of course, 
under the prior density itself, call prices were cheaper because of the lower volatility, and this effect 
seems to persist. We also see that the fatter (exponential) tails of the MED translate into higher prices 
of deeply in- or out-of-the-money call options when compared to the other two densities. 

However, as we add call prices at more strikes as constraints, the differences become smaller. In the 
third part of table [1] the prices of both call and digital options are clearly converging. 

5.2 SPX Option Prices and Heston, Schoebel-Zhu and VG Prior Densities 

In our second example, we look at call (ticker symbol SPX) and digital (BSZ) options on the Standard 
and Poor's 500 stock index [T|5]. The market data is from 18 July 2011. We consider those options 
which expire on 17 December 2011 and calibrate a Heston, Schoebel-Zhu and VG model, respectively, to 
call prices for 15 strikes 900,950, ...1550, 1600 at constant intervals of 50 using a Levenberg-Marquardt 
least-squares method ([32], [10]). The model parameters we obtained are given in table [2l 



Table 2: Model Parameters 



Parameter 


Heston 


SZ 


VG 


K 


0.8568 


1.6316 


n/a 





0.0800 


0.1731 


-0.2808 


P 


-0.8016 


-0.8031 


n/a 


a 


0.5473 


0.3249 


0.1535 


vo 


0.0421 


0.1887 


n/a 


V 


n/a 


n/a 


0.3638 



Figure [T] shows the market implied volatility skew and the volatility skews generated by the three 
models. Apart from the last two strikes at 1550 and 1600, the fit looks quite good in all three cases: 

Using formulas for the densities directly, if available, or otherwise numerical inversion (1171) . we plot 
the densities for S(T) given by these models in figure [2] The Heston and Schoebel-Zhu densities are 
almost indistinguishable from one another, whereas the VG density has a somewhat different shape with 
a slightly thinner right tail. 



Table 3: Squared Errors and Relative Entropy H(h\\p) 



Model: 


Heston 


SZ 


VG 




1.59E-04 


1.86E-04 


1.74E-04 


Relative Entropy: 


0.1311 


0.1315 


0.1374 



Table [3] shows the sum of squared errors J2i=i( (T i PX ~ <7™ odei ) 2 between the market (SPX) implied 
volatilities and the model implied volatilities. The relative entropy H{q\\p) can be seen as an alternative 
measure of fit, since by © it measures how much the prior density p needs to be deformed to obtain a 
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Figure 1: Graphs of the four volatility skews 



density q that perfectly matches the given market data. Interestingly, the Heston model fits best under 
either criterion, but the order of the Schoebel-Zhu and VG models is reversed in the two cases. 

Table 4: Digital Prices 



Strike Market (Mid) Call Spreads MED BK Heston MRED Heston SZ MRED SZ VG MRED VG 



900 





9500 





9560 





9609 





9682 





9612 





9678 





9612 





9687 


0.9607 


950 





9500 





9540 





9544 





9532 





9544 





9527 





9544 





9537 


0.9544 


1000 





9400 





9335 





9432 





9325 





9433 





9319 





9433 





9331 


0.9432 


1050 





9100 





8935 





8830 





9048 





8830 





9041 





8830 





9052 


0.8830 


1100 





8800 





8685 





8695 





8683 





8695 





8675 





8695 





8679 


0.8695 


1150 





8300 





8250 





8409 





8207 





8409 





8198 





8409 





8191 


0.8409 


1200 





7700 





7555 





7520 





7595 





7520 





7585 





7520 





7560 


0.7520 


1250 





6850 





6810 





6893 





6809 





6893 





6797 





6893 





6759 


0.6893 


1300 





5850 





5745 





5797 





5792 





5797 





5783 





5797 





5756 


0.5797 


1350 





4550 





4365 





4386 





4482 





4388 





4481 





4388 





4527 


0.4383 


1400 





3100 





2885 





2855 





2913 





2858 





2924 





2858 





3060 


0.2867 


1450 





1700 





1605 





1505 





1464 





1502 





1489 





1502 





1443 


0.1485 


1500 





0700 





0743 





0714 





0590 





0714 





0616 





0714 





0537 


0.0718 


1550 





0350 





0235 





0121 





0216 





0118 





0230 





0117 





0201 


0.0118 


1600 





0300 





0025 





0009 





0077 





0010 





0082 





0011 





0077 


0.0010 



Finally, digital prices are reported in table SJ There are noticeable differences between market prices 
(although these must be taken with a pretty big pinch of salt due to the poor liquidity and large bid-ask 
spreads), the Buchen-Kelly prices and the prices given by the three models. However, regarding the 
three relative entropy distributions obtained using the different model priors, it seems that the effect of 
the prior density on digital prices is negligible: all three MREDs basically agree with the Buchen-Kelly 
MED. 



1G 




Figure 2: Graphs of the three model densities 



6 Conclusion 



In this article we generalise the algorithm presented in [27] to the relative entropy case. The algorithm 
allows for efficient computation of a risk-neutral probability density that exactly gives European call 
option prices quoted in the market, while staying as close as possible to a given prior density under the 
criterion of relative entropy. 

It is not necessary to have an analytic expression for the prior density in question. In practice, 
several popular equity and FX models work through their characteristic functions and numerical Fourier 
inversion techniques. We pick two of these as examples, namely the Heston and the Schoebel-Zhu 
stochastic volatility models, and show how they nevertheless can be used to provide the prior density 
and incorporated into our algorithm. In other cases, analytic expressions for the density are available, 
such as for the Black-Scholes model and the Variance Gamma mode{]], and we also incorporate these 
into our analysis. 

As an application, we study the impact the choice of prior density has. In a first, purely hypothetical 
scenario, we assume that only the prices of a few options are quoted. We observe that using a prior 
density does indeed lead to significantly different option prices when compared to pricing with a pure 
log-normal density or a piece-wise exponential Buchen-Kelly density. 

In a second scenario we use option price data for S & P 500 index options for a fixed maturity traded 
on the CBOE. We calibrate three different models to these data and observe that, although the models 
generate noticeably different digital option prices, the prices obtained when using minimum relative 
entropy densities, with these models for the prior densities, agree almost perfectly. Furthermore, these 
prices are essentially the same as those given by the (non-relative) Buchen-Kelly density itself. In other 
words, in a sufficiently liquid market the effect of the prior density seems to vanish almost completely. 

1 For example, in CH — h this is easily implemented using the functions boost::math::tgamma and boost: :math::cyl_bessel_k 
in Boost [6]. 
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